A11

1 cargamos los datos

library(tidyverse)
library(patchwork)

b2018 <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/b2018.csv", header = TRUE)

b2019 <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/b2019.csv", header = TRUE)

b2020 <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/b2020.csv", header = TRUE)

bg <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/bg.csv", header = TRUE)

2 grafica de los datos

Ahora creemos 3 graficos para cada año correspondiente de nuestro conjunto de datos, cargemos algunas librerias necesarias

La siguiente forma de mapear nuestro conjunto de datos obtenido, es de la siguiente forma, considere que esta forma no involucra leer un shape file, convertir los datos a tipo geodata y todas esas particularidades, eso es lo interesante y practico, aunque no descartamos hacerlo de la forma correspondiente.

Algo particular es que si necesitamos el shape file, solo para colocar el contorno del estado de chiapas, es por ello que debemos de leerlo.

cargamos algunas librerias

library(rgdal)
library(dplyr)
library(ggplot2)
library(leaflet)      # libreria para graficar mapas interactivos
library(sf)           # manejo de informacion geografica 
library(viridis)      # paletas de colores
library(RColorBrewer) # mas paletas de colores
library(patchwork)

a continuación leemos nuestro archivo shapefile, algo que debemos identificar del shapefile es que su proyeccion geografica no coincide con la proyeccion geografica del conjunto de daos a trabajar, es por ello que se necesita realizar una transformación en la proyección de nuestro shapefile, para llevarlo a la proyección habitual de longitud y latitud que se conoce.

my_spdf <- readOGR( 
  dsn= "C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/A2", 
  layer="ENTIDAD",
  verbose=FALSE
)

# trasformamos a el sistema de coordenadas habitual

my_spdf <- spTransform(my_spdf, CRS("+proj=longlat +datum=WGS84"))

y seleccionamos el estado de chiapas

my_spdf_c <- my_spdf[my_spdf$nombre == "CHIAPAS",]

podemos observar mediante un grafico que el estado seleccionado sea el correcto.

plot(my_spdf_c, col="#f2f2f2", bg="skyblue", lwd=0.25)

Ahora si estamos listos para realizar un grafico de nuestro conjunto de datos

b2018 %>% 
  ggplot() +
  geom_point(aes(x = Longitud, y = Latitud, colour = Rain),size =3)+
  borders(my_spdf_c)+
  coord_quickmap()+
  theme_test()

Agregando un poco mas de detalle podemos obtener un grafico mas detallado para cada año y general

library(ggrepel)

b2018 %>% 
  ggplot(aes(x = Longitud, y = Latitud)) +
  borders(my_spdf_c, fill = "antiquewhite1") +
  geom_point(aes(colour = Rain), size = 4) +
  scale_color_viridis_c(option = "plasma", trans = "sqrt",
                        oob = scales::squish) +
  coord_quickmap() +
  theme_test()  +
  xlab("Longitude") + ylab("Latitude") +
  ggtitle("Chiapas 2018", subtitle = "Precipitación media por hora") +
  theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", 
        size = 0.5), panel.background = element_rect(fill = "aliceblue")) +
  geom_text_repel(data = distinct(b2018, Longitud, .keep_all = TRUE), 
                  aes(x = Longitud, y = Latitud, label = Station),
                  box.padding   = 1, point.padding = 1, segment.color = 'grey50',
                  size = 4, point.size = 4, segment.size = 1) -> p5;p5

# unique(b2018$Station)
# distinct(b2018, Longitud, .keep_all = T)
# table(b2018$Latitud)
b2019 %>% 
  ggplot() +
  borders(my_spdf_c, fill = "antiquewhite1") +
  geom_point(aes(x = Longitud, y = Latitud, color = Rain, label=Station), size =4) +
  geom_text_repel(data = distinct(b2019, Station, .keep_all = TRUE), 
                  aes(x = Longitud, y = Latitud, label = Station),
                  box.padding   = 1, point.padding = 1, segment.color = 'grey50',
                  fontface = "bold", size = 4, point.size = 4, segment.size = 1)+
  coord_quickmap() +
  theme_test()  +
  scale_color_viridis_c(option = "plasma", trans = "sqrt") +
  xlab("Longitude") + ylab("Latitude") +
  ggtitle("Chiapas 2019", subtitle = "Precipitación media por hora") +
  theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", size = 0.5), 
        panel.background = element_rect(fill = "aliceblue")) -> p6;p6

b2020 %>% 
ggplot() +
  borders(my_spdf_c, fill = "antiquewhite1") +
  geom_point(aes(x = Longitud, y = Latitud, color = Rain, label = Station), size =4) +
  geom_text_repel(data = distinct(b2020, Station, .keep_all = TRUE), 
                  aes(x = Longitud, y = Latitud, label = Station),
                  box.padding   = 1, point.padding = 1, segment.color = 'grey50',
                  fontface = "bold", size = 4, point.size = 4, segment.size = 1)+  
  coord_quickmap() +
  theme_test()  +
  scale_color_viridis_c(option = "plasma", trans = "sqrt") +
  xlab("Longitude") + ylab("Latitude") +
  ggtitle("Chiapas 2020", subtitle = "Precipitación media por hora") +
  theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", 
        size = 0.5), panel.background = element_rect(fill = "aliceblue")) -> p7;p7

bg %>% 
ggplot() +
  borders(my_spdf_c, fill = "antiquewhite1") +
  geom_point(aes(x = Longitud, y = Latitud, color = Rain, label=Station), size =4) +
  geom_text_repel(data = distinct(b2020, Station, .keep_all = TRUE), 
                  aes(x = Longitud, y = Latitud, label = Station),
                  box.padding   = 1, point.padding = 1, segment.color = 'grey50',
                  fontface = "bold", size = 4, point.size = 4, segment.size = 1)+    
  coord_quickmap() +
  theme_test()  +
  scale_color_viridis_c(option = "plasma", trans = "sqrt") +
  xlab("Longitude") + ylab("Latitude") +
  ggtitle("Chiapas 2018-2020", subtitle = "Precipitación media por hora") +
  theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", 
        size = 0.5), panel.background = element_rect(fill = "aliceblue")) -> p8;p8

(p5 + p6)/(p7 + p8)

3 Simulación del GEV

Para este caso, lo que necesitamos son ubicaciones en el estado de Chiapas, esto con el objetivo de poder simular en esas ubicaciones el modelo de valores extremos, para ello vamos utilizar la libreria geoR, esto con el fin de predecir en este caso 40 ubicaciones en la zona de estudio.

library(geoR)

bor1 <- my_spdf_c@polygons[[1]]@Polygons[[1]]@coords

sim1 <- grf(n= 40, grid = "irreg", cov.pars=c(2, 1), borders = bor1, nugget = 5, mean = 200,
            cov.model = "gaussian")

> grf: process with  1  covariance structure(s)
> grf: nugget effect is: tausq= 5 
> grf: covariance model 1 is: gaussian(sigmasq=2, phi=1)
> grf: decomposition algorithm used is:  cholesky 
> grf: End of simulation procedure. Number of realizations: 1

Podemos observar el mapa de las ubicaciones de interes simuladas, estas ubicaciones se encuentran aleatoriamente en la zona de estudio, lo unico que nos interesa en este caso es la ubicación de estos puntos simulados

par(mar=c(3,3,3,0))
points(sim1, main ="Chiapas SIM")

3.1 ajuste de la gev

Conforme a la metodologia de Davison, A.C., Padoan, S.A. and Ribatet. Lo que realizo es en los 40 sitios generados aleatoriamente, voy a obtener 50 observaciones en cada uno de ellos mediante la generación de valores de una

library(SpatialExtremes)
## Not run:
## Generate realizations from the model

n.site <- 40
n.obs <- 50

# class(sim1$coords)
coord <- cbind(lon = sim1$coords[,1], lat = sim1$coords[,2])

gp.loc   <- rgp(1, coord, "powexp", sill = 4, range = 20, smooth = 1)
gp.scale <- rgp(1, coord, "powexp", sill = 4, range = 20, smooth = 1)
gp.shape <- rgp(1, coord, "powexp", sill = 2, range = 20, smooth = 1)

# locs   <- 26 + 0.5 * coord[,"lon"] + gp.loc
# scales <- 10 + 0.2 * coord[,"lat"] + gp.scale
# shapes <- 0.15 + gp.shape

locs   <-  gp.loc
scales <-  abs(gp.scale)
shapes <-  gp.shape



data <- matrix(NA, n.obs, n.site)

for (i in 1:n.site)
  data[,i] <- rgev(n.obs, locs[i], scales[i], shapes[i])

loc.form   <- y ~ 1
scale.form <- y ~ 1
shape.form <- y ~ 1

hyper <- list()
hyper$sills <- list(loc = c(1,8), scale = c(1,1), shape = c(1,0.02))
hyper$ranges <- list(loc = c(2,20), scale = c(1,5), shape = c(1, 10))
hyper$smooths <- list(loc = c(1,1/3), scale = c(1,1/3), shape = c(1, 1/3))
hyper$betaMeans <- list(loc = 20, 
                        scale = 15, 
                        shape = 10)
hyper$betaIcov <- list(loc   = solve(diag(c(10), 1, 1)),
                       scale = solve(diag(c(10), 1, 1)),
                       shape = solve(diag(c(10), 1, 1)))

## We will use an exponential covariance function so the jump sizes for
## the shape parameter of the covariance function are null.

prop <- list(gev = c(1.2, 0.08, 0.08), ranges = c(0.7, 0.8, 0.7), smooths = c(0,0,0))
start <- list(sills   = c(4, .36, 0.009), 
              ranges  = c(24, 17, 16), 
              smooths = c(1, 1, 1), 
              beta = list(loc   = c(0.2), 
                          scale = c(0.3),
                          shape = c(20)))

# mc <- latent(data, coord, loc.form = loc.form, scale.form = scale.form,
#              shape.form = shape.form, hyper = hyper, prop = prop, start = start,
#              n = 10000, burn.in = 5000, thin = 15)
# mc

mc1 <- latent(data, coord, 
             loc.form = loc.form, 
             scale.form = scale.form,
             shape.form = shape.form, 
             hyper = hyper, 
             prop = prop, 
             start = start,
             n = 10000, 
             burn.in = 5000, 
             thin = 15)
mc1

> Effective length: 10000 
>          Burn-in: 5000 
>         Thinning: 15 
>    Effective NoP: 57.72025 
>              DIC: 10420.48 
> 
>   Regression Parameters:
>       Location Parameters:
>            lm1   
> ci.lower    7.986
> post.mean  14.230
> ci.upper   20.550
> 
>          Scale Parameters:
>            lm1  
> ci.lower   1.604
> post.mean  3.519
> ci.upper   7.058
> 
>          Shape Parameters:
>            lm1   
> ci.lower   0.2127
> post.mean  0.7418
> ci.upper   1.8975
> 
> 
>   Latent Parameters:
>       Location Parameters:
>         Covariance family: powexp 
>            sill     range    smooth 
> ci.lower     8.346   33.328    1.000
> post.mean   25.614   90.640    1.000
> ci.upper    58.684  178.641    1.000
> 
>       Scale Parameters:
>      Covariance family: powexp 
>            sill     range    smooth 
> ci.lower    0.3121   2.3281   1.0000
> post.mean   1.4136  10.4406   1.0000
> ci.upper    4.8123  26.8158   1.0000
> 
>       Shape Parameters:
>      Covariance family: powexp 
>            sill      range     smooth  
> ci.lower    0.02498   0.25977   1.00000
> post.mean   0.19374   3.59892   1.00000
> ci.upper    1.06934  19.52523   1.00000

mc2 <- latent(data, coord, 
             loc.form = loc.form, 
             scale.form = scale.form,
             shape.form = shape.form, 
             hyper = hyper, 
             prop = prop, 
             start = start,
             n = 10000, 
             burn.in = 5000, 
             thin = 15)
mc2

> Effective length: 10000 
>          Burn-in: 5000 
>         Thinning: 15 
>    Effective NoP: 57.72025 
>              DIC: 10420.48 
> 
>   Regression Parameters:
>       Location Parameters:
>            lm1   
> ci.lower    7.986
> post.mean  14.230
> ci.upper   20.550
> 
>          Scale Parameters:
>            lm1  
> ci.lower   1.604
> post.mean  3.519
> ci.upper   7.058
> 
>          Shape Parameters:
>            lm1   
> ci.lower   0.2127
> post.mean  0.7418
> ci.upper   1.8975
> 
> 
>   Latent Parameters:
>       Location Parameters:
>         Covariance family: powexp 
>            sill     range    smooth 
> ci.lower     8.346   33.328    1.000
> post.mean   25.614   90.640    1.000
> ci.upper    58.684  178.641    1.000
> 
>       Scale Parameters:
>      Covariance family: powexp 
>            sill     range    smooth 
> ci.lower    0.3121   2.3281   1.0000
> post.mean   1.4136  10.4406   1.0000
> ci.upper    4.8123  26.8158   1.0000
> 
>       Shape Parameters:
>      Covariance family: powexp 
>            sill      range     smooth  
> ci.lower    0.02498   0.25977   1.00000
> post.mean   0.19374   3.59892   1.00000
> ci.upper    1.06934  19.52523   1.00000

mc3 <- latent(data, coord, 
             loc.form = loc.form, 
             scale.form = scale.form,
             shape.form = shape.form, 
             hyper = hyper, 
             prop = prop, 
             start = start,
             n = 10000, 
             burn.in = 5000, 
             thin = 15)
mc3

> Effective length: 10000 
>          Burn-in: 5000 
>         Thinning: 15 
>    Effective NoP: 57.72025 
>              DIC: 10420.48 
> 
>   Regression Parameters:
>       Location Parameters:
>            lm1   
> ci.lower    7.986
> post.mean  14.230
> ci.upper   20.550
> 
>          Scale Parameters:
>            lm1  
> ci.lower   1.604
> post.mean  3.519
> ci.upper   7.058
> 
>          Shape Parameters:
>            lm1   
> ci.lower   0.2127
> post.mean  0.7418
> ci.upper   1.8975
> 
> 
>   Latent Parameters:
>       Location Parameters:
>         Covariance family: powexp 
>            sill     range    smooth 
> ci.lower     8.346   33.328    1.000
> post.mean   25.614   90.640    1.000
> ci.upper    58.684  178.641    1.000
> 
>       Scale Parameters:
>      Covariance family: powexp 
>            sill     range    smooth 
> ci.lower    0.3121   2.3281   1.0000
> post.mean   1.4136  10.4406   1.0000
> ci.upper    4.8123  26.8158   1.0000
> 
>       Shape Parameters:
>      Covariance family: powexp 
>            sill      range     smooth  
> ci.lower    0.02498   0.25977   1.00000
> post.mean   0.19374   3.59892   1.00000
> ci.upper    1.06934  19.52523   1.00000

dsgn.mat The design matrix.

sill = umbral

ranges = rango

smooths = forma

3.1.1 diagnosticos

Un vector de longitud 3 que devuelve el DIC, el número efectivo de parámetros eNoP y una estimación de la desviación esperada Dbar.

library(broom)
library(coda)
library(broom.mixed)
library(brms)
library(bayesplot)
library(lattice)

3.1.2 parametro de localización

library(lattice)
mc1$chain.loc <- as.data.frame(mc1$chain.loc) %>% 
  mutate(chain = 1)
mc2$chain.loc <- as.data.frame(mc2$chain.loc) %>% 
  mutate(chain = 2)
mc3$chain.loc <- as.data.frame(mc3$chain.loc) %>% 
  mutate(chain = 3)

mc.loc <-  mc1$chain.loc %>% 
  union_all(mc2$chain.loc) %>% union_all(mc3$chain.loc)

post <- posterior_samples(mc.loc, add_chain = T)
mcmc_dens_overlay(post)

post %>% 
  select(lm1, sill, range, loc1, loc2, loc3, loc4) %>% 
  mcmc_acf()

post %>% 
  select(lm1, sill, range, loc1, loc2, loc3, loc4) %>% 
  mcmc_trace()

post %>% 
  select(lm1, sill, range, loc1, loc2, loc3, loc4) %>% 
  mcmc_dens()

post %>% select(lm1) %>% mcmc_trace()

library(coda)

coda::raftery.diag(posterior_samples(post))

> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                               
>         Burn-in  Total Lower bound  Dependence
>         (M)      (N)   (Nmin)       factor (I)
>  lm1    2        3994  3746          1.07     
>  sill   3        4096  3746          1.09     
>  range  4        4673  3746          1.25     
>  smooth <NA>     <NA>  3746            NA     
>  loc1   8        10268 3746          2.74     
>  loc2   4        5080  3746          1.36     
>  loc3   8        9798  3746          2.62     
>  loc4   6        8112  3746          2.17     
>  loc5   9        9495  3746          2.53     
>  loc6   4        4752  3746          1.27     
>  loc7   3        4061  3746          1.08     
>  loc8   5        6024  3746          1.61     
>  loc9   3        4373  3746          1.17     
>  loc10  4        4913  3746          1.31     
>  loc11  9        12591 3746          3.36     
>  loc12  4        4884  3746          1.30     
>  loc13  4        4831  3746          1.29     
>  loc14  4        4831  3746          1.29     
>  loc15  4        5254  3746          1.40     
>  loc16  5        5390  3746          1.44     
>  loc17  3        4061  3746          1.08     
>  loc18  4        5038  3746          1.34     
>  loc19  6        6462  3746          1.73     
>  loc20  10       12176 3746          3.25     
>  loc21  3        4232  3746          1.13     
>  loc22  4        4996  3746          1.33     
>  loc23  8        8572  3746          2.29     
>  loc24  6        8466  3746          2.26     
>  loc25  6        7200  3746          1.92     
>  loc26  6        8462  3746          2.26     
>  loc27  8        9006  3746          2.40     
>  loc28  7        7400  3746          1.98     
>  loc29  4        4673  3746          1.25     
>  loc30  4        4673  3746          1.25     
>  loc31  6        8036  3746          2.15     
>  loc32  6        7656  3746          2.04     
>  loc33  8        11482 3746          3.07     
>  loc34  6        9066  3746          2.42     
>  loc35  15       17613 3746          4.70     
>  loc36  6        8122  3746          2.17     
>  loc37  4        4597  3746          1.23     
>  loc38  4        4831  3746          1.29     
>  loc39  4        4752  3746          1.27     
>  loc40  8        9117  3746          2.43     
>  chain  69075    69075 3746         18.40

# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.loc[,c(-4,-9)]), 
#                             as.mcmc(mc2$chain.loc[,c(-4,-9)]), 
#                             as.mcmc(mc3$chain.loc[,c(-4,-9)])))
# 
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.loc[,1]), 
#                             as.mcmc(mc2$chain.loc[,1]), 
#                             as.mcmc(mc3$chain.loc[,1])))

BMu1.mcmc<-mcmc.list(as.mcmc(mc1$chain.loc[,c(-4,-9)]), 
                     as.mcmc(mc2$chain.loc[,c(-4,-9)]),
                     as.mcmc(mc3$chain.loc[,c(-4,-9)]))

summary(BMu1.mcmc)

> 
> Iterations = 1:10000
> Thinning interval = 1 
> Number of chains = 3 
> Sample size per chain = 10000 
> 
> 1. Empirical mean and standard deviation for each variable,
>    plus standard error of the mean:
> 
>          Mean      SD  Naive SE Time-series SE
> lm1   14.2297  3.2277 0.0186353       0.020648
> sill  25.6144 12.8728 0.0743212       0.093081
> range 90.6396 37.5718 0.2169210       0.268302
> loc1   1.6649  0.1669 0.0009634       0.001768
> loc2   1.1084  0.2365 0.0013654       0.002141
> loc3   2.6815  0.1991 0.0011494       0.002124
> loc4   1.4826  0.1550 0.0008947       0.001900
> loc6   1.5297  0.2102 0.0012134       0.002090
> loc7   2.3124  0.3217 0.0018574       0.002499
> loc8   1.1358  0.2054 0.0011857       0.001996
> loc9   2.2184  0.1945 0.0011229       0.001706
> loc10  2.0980  0.1892 0.0010926       0.001660
> loc11  1.8553  0.1551 0.0008955       0.001669
> loc12  2.8993  0.2871 0.0016574       0.002579
> loc13  1.3152  0.1928 0.0011133       0.001533
> loc14  2.1239  0.2245 0.0012962       0.001910
> loc15  0.9606  0.1953 0.0011275       0.001722
> loc16  0.8013  0.2432 0.0014039       0.001933
> loc17  2.0278  0.1857 0.0010719       0.001463
> loc18  1.5992  0.2115 0.0012212       0.002068
> loc19  1.7372  0.1574 0.0009088       0.001827
> loc20  1.9681  0.1521 0.0008783       0.002631
> loc21  1.8076  0.2176 0.0012565       0.001625
> loc22  2.7312  0.2012 0.0011615       0.001811
> loc23  2.4389  0.2191 0.0012649       0.002205
> loc24  1.4109  0.2311 0.0013343       0.002245
> loc25  1.7159  0.1478 0.0008534       0.001889
> loc26  1.4967  0.2041 0.0011787       0.001733
> loc27  1.4447  0.1708 0.0009860       0.001962
> loc28  1.7872  0.1631 0.0009415       0.002100
> loc29  1.2776  0.3351 0.0019347       0.002832
> loc30  2.6040  0.2886 0.0016663       0.002571
> loc31  1.3298  0.1946 0.0011238       0.001751
> loc32  1.6352  0.1847 0.0010664       0.002063
> loc33  1.7319  0.1546 0.0008929       0.001864
> loc34  1.6563  0.1839 0.0010620       0.001982
> loc35  1.9140  0.1555 0.0008977       0.002602
> loc36  2.9123  0.1927 0.0011128       0.001833
> loc37  2.7145  0.2101 0.0012130       0.001996
> loc38  1.4030  0.1941 0.0011207       0.001710
> loc39  0.7753  0.2160 0.0012473       0.001728
> loc40  1.6935  0.1554 0.0008973       0.001975
> chain  2.0000  0.8165 0.0047141       0.000000
> 
> 2. Quantiles for each variable:
> 
>          2.5%     25%     50%      75%   97.5%
> lm1    7.9863 12.0235 14.1928  16.4447  20.550
> sill   8.3460 16.5638 22.9614  31.7648  58.684
> range 33.3278 63.4518 85.1777 111.7662 178.641
> loc1   1.3388  1.5532  1.6633   1.7754   2.005
> loc2   0.6341  0.9545  1.1121   1.2688   1.564
> loc3   2.3052  2.5479  2.6767   2.8095   3.082
> loc4   1.1857  1.3775  1.4792   1.5869   1.794
> loc6   1.1307  1.3868  1.5268   1.6691   1.955
> loc7   1.7000  2.0965  2.3065   2.5208   2.961
> loc8   0.7305  0.9979  1.1367   1.2741   1.535
> loc9   1.8506  2.0848  2.2127   2.3483   2.609
> loc10  1.7318  1.9717  2.0970   2.2243   2.480
> loc11  1.5622  1.7490  1.8529   1.9575   2.167
> loc12  2.3618  2.7045  2.8911   3.0802   3.503
> loc13  0.9393  1.1859  1.3147   1.4411   1.705
> loc14  1.6874  1.9732  2.1214   2.2720   2.573
> loc15  0.5750  0.8288  0.9608   1.0926   1.343
> loc16  0.3258  0.6391  0.7956   0.9628   1.285
> loc17  1.6709  1.9015  2.0237   2.1529   2.392
> loc18  1.1934  1.4570  1.5992   1.7395   2.024
> loc19  1.4373  1.6305  1.7339   1.8427   2.052
> loc20  1.6848  1.8640  1.9645   2.0684   2.279
> loc21  1.3817  1.6618  1.8072   1.9536   2.233
> loc22  2.3526  2.5925  2.7260   2.8648   3.140
> loc23  2.0164  2.2882  2.4370   2.5840   2.879
> loc24  0.9653  1.2586  1.4071   1.5629   1.865
> loc25  1.4339  1.6136  1.7136   1.8143   2.014
> loc26  1.1085  1.3571  1.4948   1.6317   1.913
> loc27  1.1114  1.3305  1.4411   1.5586   1.782
> loc28  1.4788  1.6756  1.7821   1.8958   2.125
> loc29  0.6158  1.0549  1.2772   1.5007   1.942
> loc30  2.0468  2.4090  2.5997   2.7892   3.180
> loc31  0.9537  1.1987  1.3286   1.4592   1.719
> loc32  1.2736  1.5121  1.6367   1.7608   2.004
> loc33  1.4307  1.6275  1.7295   1.8322   2.040
> loc34  1.2974  1.5318  1.6559   1.7798   2.015
> loc35  1.6182  1.8086  1.9090   2.0170   2.224
> loc36  2.5401  2.7812  2.9073   3.0403   3.304
> loc37  2.3147  2.5724  2.7103   2.8543   3.135
> loc38  1.0310  1.2713  1.4035   1.5297   1.791
> loc39  0.3605  0.6279  0.7756   0.9167   1.213
> loc40  1.3970  1.5887  1.6923   1.7956   2.007
> chain  1.0000  1.0000  2.0000   3.0000   3.000

xyplot(BMu1.mcmc)

densityplot(BMu1.mcmc)                             #Densidades

plot(BMu1.mcmc)

# layout(matrix(1:12, 3,4));  

traceplot(BMu1.mcmc)

#Grafica de autocorrelacion
autocorr.plot(BMu1.mcmc, auto.layout = TRUE, ask =F)

coda::acfplot(BMu1.mcmc)

# gelman.diag(BMu1.mcmc)
# 
# gelman.plot(BMu1.mcmc)

geweke.diag(BMu1.mcmc)

> [[1]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>      lm1     sill    range     loc1     loc2     loc3     loc4     loc6 
> -2.02360 -0.13114 -0.57292 -0.13784  1.40118 -0.25933  0.79576  0.38023 
>     loc7     loc8     loc9    loc10    loc11    loc12    loc13    loc14 
> -0.65418 -0.81671  0.76711  1.53365 -0.43782 -0.78151  0.43533  0.83240 
>    loc15    loc16    loc17    loc18    loc19    loc20    loc21    loc22 
>  0.21704 -0.02015  2.09454  0.25178 -1.06965 -0.58544  1.15215 -2.70446 
>    loc23    loc24    loc25    loc26    loc27    loc28    loc29    loc30 
> -1.37479 -0.25433 -0.39967 -0.52450  1.25959 -1.75285 -0.62833 -1.14631 
>    loc31    loc32    loc33    loc34    loc35    loc36    loc37    loc38 
> -0.50914  0.96571 -0.59595  0.91762 -0.59011  0.74688  0.20454 -0.78983 
>    loc39    loc40    chain 
> -1.10704 -0.34407      NaN 
> 
> 
> [[2]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>      lm1     sill    range     loc1     loc2     loc3     loc4     loc6 
> -2.02360 -0.13114 -0.57292 -0.13784  1.40118 -0.25933  0.79576  0.38023 
>     loc7     loc8     loc9    loc10    loc11    loc12    loc13    loc14 
> -0.65418 -0.81671  0.76711  1.53365 -0.43782 -0.78151  0.43533  0.83240 
>    loc15    loc16    loc17    loc18    loc19    loc20    loc21    loc22 
>  0.21704 -0.02015  2.09454  0.25178 -1.06965 -0.58544  1.15215 -2.70446 
>    loc23    loc24    loc25    loc26    loc27    loc28    loc29    loc30 
> -1.37479 -0.25433 -0.39967 -0.52450  1.25959 -1.75285 -0.62833 -1.14631 
>    loc31    loc32    loc33    loc34    loc35    loc36    loc37    loc38 
> -0.50914  0.96571 -0.59595  0.91762 -0.59011  0.74688  0.20454 -0.78983 
>    loc39    loc40    chain 
> -1.10704 -0.34407      NaN 
> 
> 
> [[3]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>      lm1     sill    range     loc1     loc2     loc3     loc4     loc6 
> -2.02360 -0.13114 -0.57292 -0.13784  1.40118 -0.25933  0.79576  0.38023 
>     loc7     loc8     loc9    loc10    loc11    loc12    loc13    loc14 
> -0.65418 -0.81671  0.76711  1.53365 -0.43782 -0.78151  0.43533  0.83240 
>    loc15    loc16    loc17    loc18    loc19    loc20    loc21    loc22 
>  0.21704 -0.02015  2.09454  0.25178 -1.06965 -0.58544  1.15215 -2.70446 
>    loc23    loc24    loc25    loc26    loc27    loc28    loc29    loc30 
> -1.37479 -0.25433 -0.39967 -0.52450  1.25959 -1.75285 -0.62833 -1.14631 
>    loc31    loc32    loc33    loc34    loc35    loc36    loc37    loc38 
> -0.50914  0.96571 -0.59595  0.91762 -0.59011  0.74688  0.20454 -0.78983 
>    loc39    loc40    chain 
> -1.10704 -0.34407      NaN

#geweke.plot(BMu1.mcmc)
raftery.diag(BMu1.mcmc)

> [[1]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                              
>        Burn-in  Total Lower bound  Dependence
>        (M)      (N)   (Nmin)       factor (I)
>  lm1   2        3994  3746         1.07      
>  sill  3        4098  3746         1.09      
>  range 4        4674  3746         1.25      
>  loc1  8        10270 3746         2.74      
>  loc2  4        5080  3746         1.36      
>  loc3  6        6349  3746         1.69      
>  loc4  6        6877  3746         1.84      
>  loc6  4        4752  3746         1.27      
>  loc7  3        4061  3746         1.08      
>  loc8  5        6025  3746         1.61      
>  loc9  3        4374  3746         1.17      
>  loc10 4        4913  3746         1.31      
>  loc11 4        4913  3746         1.31      
>  loc12 4        4884  3746         1.30      
>  loc13 4        4832  3746         1.29      
>  loc14 4        4832  3746         1.29      
>  loc15 4        5254  3746         1.40      
>  loc16 5        5390  3746         1.44      
>  loc17 3        4061  3746         1.08      
>  loc18 4        5038  3746         1.34      
>  loc19 6        6462  3746         1.73      
>  loc20 9        9691  3746         2.59      
>  loc21 3        4232  3746         1.13      
>  loc22 4        4996  3746         1.33      
>  loc23 5        6247  3746         1.67      
>  loc24 4        5166  3746         1.38      
>  loc25 6        7200  3746         1.92      
>  loc26 3        4520  3746         1.21      
>  loc27 6        6816  3746         1.82      
>  loc28 7        7401  3746         1.98      
>  loc29 4        4674  3746         1.25      
>  loc30 4        4674  3746         1.25      
>  loc31 4        5080  3746         1.36      
>  loc32 5        5820  3746         1.55      
>  loc33 7        7465  3746         1.99      
>  loc34 4        5254  3746         1.40      
>  loc35 10       11232 3746         3.00      
>  loc36 3        4520  3746         1.21      
>  loc37 4        4597  3746         1.23      
>  loc38 4        4832  3746         1.29      
>  loc39 4        4752  3746         1.27      
>  loc40 8        9118  3746         2.43      
>  chain <NA>     <NA>  3746           NA      
> 
> 
> [[2]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                              
>        Burn-in  Total Lower bound  Dependence
>        (M)      (N)   (Nmin)       factor (I)
>  lm1   2        3994  3746         1.07      
>  sill  3        4098  3746         1.09      
>  range 4        4674  3746         1.25      
>  loc1  8        10270 3746         2.74      
>  loc2  4        5080  3746         1.36      
>  loc3  6        6349  3746         1.69      
>  loc4  6        6877  3746         1.84      
>  loc6  4        4752  3746         1.27      
>  loc7  3        4061  3746         1.08      
>  loc8  5        6025  3746         1.61      
>  loc9  3        4374  3746         1.17      
>  loc10 4        4913  3746         1.31      
>  loc11 4        4913  3746         1.31      
>  loc12 4        4884  3746         1.30      
>  loc13 4        4832  3746         1.29      
>  loc14 4        4832  3746         1.29      
>  loc15 4        5254  3746         1.40      
>  loc16 5        5390  3746         1.44      
>  loc17 3        4061  3746         1.08      
>  loc18 4        5038  3746         1.34      
>  loc19 6        6462  3746         1.73      
>  loc20 9        9691  3746         2.59      
>  loc21 3        4232  3746         1.13      
>  loc22 4        4996  3746         1.33      
>  loc23 5        6247  3746         1.67      
>  loc24 4        5166  3746         1.38      
>  loc25 6        7200  3746         1.92      
>  loc26 3        4520  3746         1.21      
>  loc27 6        6816  3746         1.82      
>  loc28 7        7401  3746         1.98      
>  loc29 4        4674  3746         1.25      
>  loc30 4        4674  3746         1.25      
>  loc31 4        5080  3746         1.36      
>  loc32 5        5820  3746         1.55      
>  loc33 7        7465  3746         1.99      
>  loc34 4        5254  3746         1.40      
>  loc35 10       11232 3746         3.00      
>  loc36 3        4520  3746         1.21      
>  loc37 4        4597  3746         1.23      
>  loc38 4        4832  3746         1.29      
>  loc39 4        4752  3746         1.27      
>  loc40 8        9118  3746         2.43      
>  chain <NA>     <NA>  3746           NA      
> 
> 
> [[3]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                              
>        Burn-in  Total Lower bound  Dependence
>        (M)      (N)   (Nmin)       factor (I)
>  lm1   2        3994  3746         1.07      
>  sill  3        4098  3746         1.09      
>  range 4        4674  3746         1.25      
>  loc1  8        10270 3746         2.74      
>  loc2  4        5080  3746         1.36      
>  loc3  6        6349  3746         1.69      
>  loc4  6        6877  3746         1.84      
>  loc6  4        4752  3746         1.27      
>  loc7  3        4061  3746         1.08      
>  loc8  5        6025  3746         1.61      
>  loc9  3        4374  3746         1.17      
>  loc10 4        4913  3746         1.31      
>  loc11 4        4913  3746         1.31      
>  loc12 4        4884  3746         1.30      
>  loc13 4        4832  3746         1.29      
>  loc14 4        4832  3746         1.29      
>  loc15 4        5254  3746         1.40      
>  loc16 5        5390  3746         1.44      
>  loc17 3        4061  3746         1.08      
>  loc18 4        5038  3746         1.34      
>  loc19 6        6462  3746         1.73      
>  loc20 9        9691  3746         2.59      
>  loc21 3        4232  3746         1.13      
>  loc22 4        4996  3746         1.33      
>  loc23 5        6247  3746         1.67      
>  loc24 4        5166  3746         1.38      
>  loc25 6        7200  3746         1.92      
>  loc26 3        4520  3746         1.21      
>  loc27 6        6816  3746         1.82      
>  loc28 7        7401  3746         1.98      
>  loc29 4        4674  3746         1.25      
>  loc30 4        4674  3746         1.25      
>  loc31 4        5080  3746         1.36      
>  loc32 5        5820  3746         1.55      
>  loc33 7        7465  3746         1.99      
>  loc34 4        5254  3746         1.40      
>  loc35 10       11232 3746         3.00      
>  loc36 3        4520  3746         1.21      
>  loc37 4        4597  3746         1.23      
>  loc38 4        4832  3746         1.29      
>  loc39 4        4752  3746         1.27      
>  loc40 8        9118  3746         2.43      
>  chain <NA>     <NA>  3746           NA

# heidel.diag(BMu1.mcmc)

3.1.3 parametro de escala

mc1$chain.scale <- as.data.frame(mc1$chain.scale) %>% 
  mutate(chain = 1)
mc2$chain.scale <- as.data.frame(mc2$chain.scale) %>% 
  mutate(chain = 2)
mc3$chain.scale <- as.data.frame(mc3$chain.scale) %>% 
  mutate(chain = 3)

mc.loc <-  mc1$chain.scale %>% 
  union_all(mc2$chain.scale) %>% union_all(mc3$chain.scale)

post <- posterior_samples(mc.loc, add_chain = T)
mcmc_dens_overlay(post)

post %>% 
  select(lm1, sill, range, scale1, scale2, scale3, scale4) %>% 
  mcmc_acf()

post %>% 
  select(lm1, sill, range, scale1, scale2, scale3, scale4) %>% 
  mcmc_trace()

post %>% 
  select(lm1, sill, range, scale1, scale2, scale3, scale4) %>% 
  mcmc_dens()

post %>% select(lm1) %>% mcmc_trace()

library(coda)

coda::raftery.diag(posterior_samples(post))

> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     2        3740  3746          0.998    
>  sill    3        4232  3746          1.130    
>  range   8        9716  3746          2.590    
>  smooth  <NA>     <NA>  3746             NA    
>  scale1  6        8308  3746          2.220    
>  scale2  3        4095  3746          1.090    
>  scale3  6        9070  3746          2.420    
>  scale4  6        8482  3746          2.260    
>  scale5  8        9888  3746          2.640    
>  scale6  6        11691 3746          3.120    
>  scale7  6        7762  3746          2.070    
>  scale8  6        7760  3746          2.070    
>  scale9  3        4520  3746          1.210    
>  scale10 6        8380  3746          2.240    
>  scale11 6        9414  3746          2.510    
>  scale12 4        5038  3746          1.340    
>  scale13 4        7048  3746          1.880    
>  scale14 3        4128  3746          1.100    
>  scale15 3        4558  3746          1.220    
>  scale16 6        8826  3746          2.360    
>  scale17 4        8092  3746          2.160    
>  scale18 8        9270  3746          2.470    
>  scale19 6        8574  3746          2.290    
>  scale20 9        13014 3746          3.470    
>  scale21 3        4337  3746          1.160    
>  scale22 3        4302  3746          1.150    
>  scale23 4        4635  3746          1.240    
>  scale24 6        8820  3746          2.350    
>  scale25 6        8826  3746          2.360    
>  scale26 4        4752  3746          1.270    
>  scale27 6        7818  3746          2.090    
>  scale28 9        12801 3746          3.420    
>  scale29 6        8190  3746          2.190    
>  scale30 4        7702  3746          2.060    
>  scale31 3        4560  3746          1.220    
>  scale32 6        7708  3746          2.060    
>  scale33 6        9016  3746          2.410    
>  scale34 6        8978  3746          2.400    
>  scale35 12       14139 3746          3.770    
>  scale36 6        7626  3746          2.040    
>  scale37 6        8278  3746          2.210    
>  scale38 6        7260  3746          1.940    
>  scale39 3        4373  3746          1.170    
>  scale40 9        13674 3746          3.650    
>  chain   69075    69075 3746         18.400

# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.scale[,c(-4,-9)]), 
#                             as.mcmc(mc2$chain.scale[,c(-4,-9)]), 
#                             as.mcmc(mc3$chain.scale[,c(-4,-9)])))
# 
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.scale[,1]), 
#                             as.mcmc(mc2$chain.scale[,1]), 
#                             as.mcmc(mc3$chain.scale[,1])))

BMu1.mcmc<-mcmc.list(as.mcmc(mc1$chain.scale[,c(-4,-9)]), 
                     as.mcmc(mc2$chain.scale[,c(-4,-9)]),
                     as.mcmc(mc3$chain.scale[,c(-4,-9)]))

summary(BMu1.mcmc)

> 
> Iterations = 1:10000
> Thinning interval = 1 
> Number of chains = 3 
> Sample size per chain = 10000 
> 
> 1. Empirical mean and standard deviation for each variable,
>    plus standard error of the mean:
> 
>           Mean     SD  Naive SE Time-series SE
> lm1      3.519 1.3910 0.0080307       0.010806
> sill     1.414 1.1966 0.0069083       0.011339
> range   10.441 6.4762 0.0373902       0.063922
> scale1   1.809 0.1598 0.0009226       0.001661
> scale2   2.207 0.1706 0.0009850       0.001530
> scale3   2.088 0.1744 0.0010071       0.001749
> scale4   1.780 0.1572 0.0009075       0.001829
> scale6   2.330 0.2129 0.0012291       0.002264
> scale7   2.901 0.2497 0.0014415       0.002190
> scale8   2.107 0.1737 0.0010027       0.001542
> scale9   2.060 0.1753 0.0010122       0.001552
> scale10  2.136 0.1762 0.0010176       0.001538
> scale11  1.610 0.1552 0.0008963       0.001704
> scale12  2.788 0.3168 0.0018289       0.002945
> scale13  1.932 0.1785 0.0010305       0.001501
> scale14  2.227 0.1890 0.0010912       0.001587
> scale15  1.870 0.1880 0.0010855       0.001704
> scale16  2.003 0.2167 0.0012511       0.001773
> scale17  1.875 0.1783 0.0010294       0.001582
> scale18  2.228 0.2141 0.0012361       0.002856
> scale19  1.689 0.1541 0.0008896       0.001689
> scale20  1.597 0.1923 0.0011104       0.003180
> scale21  2.187 0.2067 0.0011936       0.001578
> scale22  1.791 0.1940 0.0011203       0.001784
> scale23  1.895 0.2030 0.0011722       0.002175
> scale24  2.380 0.2196 0.0012676       0.002243
> scale25  1.549 0.1516 0.0008752       0.001925
> scale26  2.051 0.1780 0.0010275       0.001535
> scale27  1.868 0.1607 0.0009279       0.001863
> scale28  1.538 0.1741 0.0010051       0.002257
> scale29  2.827 0.2687 0.0015512       0.002566
> scale30  2.618 0.2675 0.0015446       0.002603
> scale31  2.083 0.1739 0.0010040       0.001541
> scale32  2.167 0.1487 0.0008585       0.001525
> scale33  1.542 0.1595 0.0009211       0.001994
> scale34  2.125 0.1563 0.0009027       0.001479
> scale35  1.596 0.1905 0.0011001       0.003125
> scale36  1.871 0.1992 0.0011501       0.001985
> scale37  2.101 0.2121 0.0012246       0.002012
> scale38  2.037 0.1885 0.0010886       0.001756
> scale39  1.746 0.1884 0.0010876       0.001579
> scale40  1.526 0.1812 0.0010464       0.002234
> chain    2.000 0.8165 0.0047141       0.000000
> 
> 2. Quantiles for each variable:
> 
>           2.5%    25%   50%    75%  97.5%
> lm1     1.6039 2.5820 3.227  4.120  7.058
> sill    0.3121 0.6665 1.057  1.711  4.812
> range   2.3281 5.7278 9.019 13.476 26.816
> scale1  1.5099 1.7028 1.806  1.911  2.133
> scale2  1.8930 2.0905 2.198  2.315  2.568
> scale3  1.7694 1.9665 2.081  2.205  2.450
> scale4  1.4845 1.6729 1.778  1.881  2.097
> scale6  1.9246 2.1867 2.323  2.465  2.761
> scale7  2.4643 2.7240 2.887  3.055  3.448
> scale8  1.7847 1.9877 2.104  2.220  2.459
> scale9  1.7218 1.9424 2.059  2.177  2.406
> scale10 1.7938 2.0180 2.132  2.252  2.497
> scale11 1.3121 1.5065 1.608  1.712  1.918
> scale12 2.2050 2.5693 2.769  2.989  3.448
> scale13 1.5807 1.8126 1.933  2.052  2.281
> scale14 1.8770 2.0964 2.220  2.352  2.617
> scale15 1.5099 1.7473 1.865  1.994  2.248
> scale16 1.6073 1.8528 1.996  2.145  2.450
> scale17 1.5136 1.7569 1.878  1.992  2.225
> scale18 1.8129 2.0853 2.223  2.368  2.653
> scale19 1.3903 1.5835 1.687  1.789  1.998
> scale20 1.2398 1.4654 1.593  1.721  2.000
> scale21 1.7864 2.0488 2.185  2.325  2.597
> scale22 1.4287 1.6574 1.785  1.916  2.183
> scale23 1.5234 1.7539 1.883  2.027  2.315
> scale24 1.9740 2.2292 2.371  2.522  2.826
> scale25 1.2591 1.4452 1.545  1.649  1.856
> scale26 1.7034 1.9331 2.047  2.167  2.405
> scale27 1.5682 1.7597 1.863  1.971  2.194
> scale28 1.2143 1.4180 1.531  1.650  1.900
> scale29 2.3530 2.6389 2.807  2.995  3.407
> scale30 2.1277 2.4370 2.604  2.787  3.188
> scale31 1.7489 1.9637 2.082  2.200  2.430
> scale32 1.8911 2.0638 2.164  2.262  2.466
> scale33 1.2348 1.4332 1.538  1.647  1.868
> scale34 1.8289 2.0176 2.122  2.228  2.441
> scale35 1.2396 1.4673 1.590  1.717  1.995
> scale36 1.4885 1.7351 1.870  2.004  2.267
> scale37 1.6931 1.9580 2.098  2.241  2.525
> scale38 1.6741 1.9092 2.035  2.162  2.414
> scale39 1.4048 1.6172 1.737  1.865  2.139
> scale40 1.1760 1.4047 1.523  1.643  1.891
> chain   1.0000 1.0000 2.000  3.000  3.000

xyplot(BMu1.mcmc)

densityplot(BMu1.mcmc)                             #Densidades

plot(BMu1.mcmc)

# layout(matrix(1:12, 3,4));  

traceplot(BMu1.mcmc)

#Grafica de autocorrelacion
autocorr.plot(BMu1.mcmc, auto.layout = TRUE, ask =F)

coda::acfplot(BMu1.mcmc)

 # gelman.diag(BMu1.mcmc)
# 
# gelman.plot(BMu1.mcmc)

geweke.diag(BMu1.mcmc)

> [[1]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>      lm1     sill    range   scale1   scale2   scale3   scale4   scale6 
> -0.08032 -0.02250  1.17062 -0.50661  1.31192  0.31638  0.87023 -0.04461 
>   scale7   scale8   scale9  scale10  scale11  scale12  scale13  scale14 
> -0.73503 -1.30912  1.91126  2.23194  0.36890 -0.70658  1.04083  1.01627 
>  scale15  scale16  scale17  scale18  scale19  scale20  scale21  scale22 
>  0.68221  0.73472  1.83715  0.09793 -0.28420 -0.44112  1.22216 -3.26133 
>  scale23  scale24  scale25  scale26  scale27  scale28  scale29  scale30 
> -2.80187  0.04558 -0.13999 -0.41964  1.43725 -1.79599 -1.00720 -0.25150 
>  scale31  scale32  scale33  scale34  scale35  scale36  scale37  scale38 
> -0.66973  1.55224 -0.21793  1.63241 -0.59266  0.97689  0.18939 -0.53277 
>  scale39  scale40    chain 
>  0.69119 -0.37436      NaN 
> 
> 
> [[2]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>      lm1     sill    range   scale1   scale2   scale3   scale4   scale6 
> -0.08032 -0.02250  1.17062 -0.50661  1.31192  0.31638  0.87023 -0.04461 
>   scale7   scale8   scale9  scale10  scale11  scale12  scale13  scale14 
> -0.73503 -1.30912  1.91126  2.23194  0.36890 -0.70658  1.04083  1.01627 
>  scale15  scale16  scale17  scale18  scale19  scale20  scale21  scale22 
>  0.68221  0.73472  1.83715  0.09793 -0.28420 -0.44112  1.22216 -3.26133 
>  scale23  scale24  scale25  scale26  scale27  scale28  scale29  scale30 
> -2.80187  0.04558 -0.13999 -0.41964  1.43725 -1.79599 -1.00720 -0.25150 
>  scale31  scale32  scale33  scale34  scale35  scale36  scale37  scale38 
> -0.66973  1.55224 -0.21793  1.63241 -0.59266  0.97689  0.18939 -0.53277 
>  scale39  scale40    chain 
>  0.69119 -0.37436      NaN 
> 
> 
> [[3]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>      lm1     sill    range   scale1   scale2   scale3   scale4   scale6 
> -0.08032 -0.02250  1.17062 -0.50661  1.31192  0.31638  0.87023 -0.04461 
>   scale7   scale8   scale9  scale10  scale11  scale12  scale13  scale14 
> -0.73503 -1.30912  1.91126  2.23194  0.36890 -0.70658  1.04083  1.01627 
>  scale15  scale16  scale17  scale18  scale19  scale20  scale21  scale22 
>  0.68221  0.73472  1.83715  0.09793 -0.28420 -0.44112  1.22216 -3.26133 
>  scale23  scale24  scale25  scale26  scale27  scale28  scale29  scale30 
> -2.80187  0.04558 -0.13999 -0.41964  1.43725 -1.79599 -1.00720 -0.25150 
>  scale31  scale32  scale33  scale34  scale35  scale36  scale37  scale38 
> -0.66973  1.55224 -0.21793  1.63241 -0.59266  0.97689  0.18939 -0.53277 
>  scale39  scale40    chain 
>  0.69119 -0.37436      NaN

#geweke.plot(BMu1.mcmc)
raftery.diag(BMu1.mcmc)

> [[1]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     2        3741  3746         0.999     
>  sill    3        4232  3746         1.130     
>  range   5        5721  3746         1.530     
>  scale1  4        4752  3746         1.270     
>  scale2  3        4095  3746         1.090     
>  scale3  4        4635  3746         1.240     
>  scale4  5        5390  3746         1.440     
>  scale6  4        4674  3746         1.250     
>  scale7  3        4320  3746         1.150     
>  scale8  4        4635  3746         1.240     
>  scale9  3        4520  3746         1.210     
>  scale10 6        8382  3746         2.240     
>  scale11 6        9414  3746         2.510     
>  scale12 4        5038  3746         1.340     
>  scale13 4        4635  3746         1.240     
>  scale14 3        4129  3746         1.100     
>  scale15 3        4558  3746         1.220     
>  scale16 3        4520  3746         1.210     
>  scale17 3        4410  3746         1.180     
>  scale18 8        9272  3746         2.480     
>  scale19 6        8574  3746         2.290     
>  scale20 6        9466  3746         2.530     
>  scale21 3        4338  3746         1.160     
>  scale22 3        4302  3746         1.150     
>  scale23 4        4635  3746         1.240     
>  scale24 4        4635  3746         1.240     
>  scale25 6        8828  3746         2.360     
>  scale26 4        4752  3746         1.270     
>  scale27 6        7818  3746         2.090     
>  scale28 8        8872  3746         2.370     
>  scale29 3        4302  3746         1.150     
>  scale30 3        4483  3746         1.200     
>  scale31 3        4564  3746         1.220     
>  scale32 3        4410  3746         1.180     
>  scale33 6        9018  3746         2.410     
>  scale34 3        4410  3746         1.180     
>  scale35 10       12532 3746         3.350     
>  scale36 4        4832  3746         1.290     
>  scale37 4        4872  3746         1.300     
>  scale38 4        4674  3746         1.250     
>  scale39 3        4374  3746         1.170     
>  scale40 6        10278 3746         2.740     
>  chain   <NA>     <NA>  3746            NA     
> 
> 
> [[2]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     2        3741  3746         0.999     
>  sill    3        4232  3746         1.130     
>  range   5        5721  3746         1.530     
>  scale1  4        4752  3746         1.270     
>  scale2  3        4095  3746         1.090     
>  scale3  4        4635  3746         1.240     
>  scale4  5        5390  3746         1.440     
>  scale6  4        4674  3746         1.250     
>  scale7  3        4320  3746         1.150     
>  scale8  4        4635  3746         1.240     
>  scale9  3        4520  3746         1.210     
>  scale10 6        8382  3746         2.240     
>  scale11 6        9414  3746         2.510     
>  scale12 4        5038  3746         1.340     
>  scale13 4        4635  3746         1.240     
>  scale14 3        4129  3746         1.100     
>  scale15 3        4558  3746         1.220     
>  scale16 3        4520  3746         1.210     
>  scale17 3        4410  3746         1.180     
>  scale18 8        9272  3746         2.480     
>  scale19 6        8574  3746         2.290     
>  scale20 6        9466  3746         2.530     
>  scale21 3        4338  3746         1.160     
>  scale22 3        4302  3746         1.150     
>  scale23 4        4635  3746         1.240     
>  scale24 4        4635  3746         1.240     
>  scale25 6        8828  3746         2.360     
>  scale26 4        4752  3746         1.270     
>  scale27 6        7818  3746         2.090     
>  scale28 8        8872  3746         2.370     
>  scale29 3        4302  3746         1.150     
>  scale30 3        4483  3746         1.200     
>  scale31 3        4564  3746         1.220     
>  scale32 3        4410  3746         1.180     
>  scale33 6        9018  3746         2.410     
>  scale34 3        4410  3746         1.180     
>  scale35 10       12532 3746         3.350     
>  scale36 4        4832  3746         1.290     
>  scale37 4        4872  3746         1.300     
>  scale38 4        4674  3746         1.250     
>  scale39 3        4374  3746         1.170     
>  scale40 6        10278 3746         2.740     
>  chain   <NA>     <NA>  3746            NA     
> 
> 
> [[3]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     2        3741  3746         0.999     
>  sill    3        4232  3746         1.130     
>  range   5        5721  3746         1.530     
>  scale1  4        4752  3746         1.270     
>  scale2  3        4095  3746         1.090     
>  scale3  4        4635  3746         1.240     
>  scale4  5        5390  3746         1.440     
>  scale6  4        4674  3746         1.250     
>  scale7  3        4320  3746         1.150     
>  scale8  4        4635  3746         1.240     
>  scale9  3        4520  3746         1.210     
>  scale10 6        8382  3746         2.240     
>  scale11 6        9414  3746         2.510     
>  scale12 4        5038  3746         1.340     
>  scale13 4        4635  3746         1.240     
>  scale14 3        4129  3746         1.100     
>  scale15 3        4558  3746         1.220     
>  scale16 3        4520  3746         1.210     
>  scale17 3        4410  3746         1.180     
>  scale18 8        9272  3746         2.480     
>  scale19 6        8574  3746         2.290     
>  scale20 6        9466  3746         2.530     
>  scale21 3        4338  3746         1.160     
>  scale22 3        4302  3746         1.150     
>  scale23 4        4635  3746         1.240     
>  scale24 4        4635  3746         1.240     
>  scale25 6        8828  3746         2.360     
>  scale26 4        4752  3746         1.270     
>  scale27 6        7818  3746         2.090     
>  scale28 8        8872  3746         2.370     
>  scale29 3        4302  3746         1.150     
>  scale30 3        4483  3746         1.200     
>  scale31 3        4564  3746         1.220     
>  scale32 3        4410  3746         1.180     
>  scale33 6        9018  3746         2.410     
>  scale34 3        4410  3746         1.180     
>  scale35 10       12532 3746         3.350     
>  scale36 4        4832  3746         1.290     
>  scale37 4        4872  3746         1.300     
>  scale38 4        4674  3746         1.250     
>  scale39 3        4374  3746         1.170     
>  scale40 6        10278 3746         2.740     
>  chain   <NA>     <NA>  3746            NA

# heidel.diag(BMu1.mcmc)

3.1.4 parametro de forma

mc1$chain.shape <- as.data.frame(mc1$chain.shape) %>% 
  mutate(chain = 1)
mc2$chain.shape <- as.data.frame(mc2$chain.shape) %>% 
  mutate(chain = 2)
mc3$chain.shape <- as.data.frame(mc3$chain.shape) %>% 
  mutate(chain = 3)

mc.loc <-  mc1$chain.shape %>% 
  union_all(mc2$chain.shape) %>% union_all(mc3$chain.shape)

post <- posterior_samples(mc.loc, add_chain = T)
mcmc_dens_overlay(post)

post %>% 
  select(lm1, sill, range, shape1, shape2, shape3, shape4) %>% 
  mcmc_acf()

post %>% 
  select(lm1, sill, range, shape1, shape2, shape3, shape4) %>% 
  mcmc_trace()

post %>% 
  select(lm1, sill, range, shape1, shape2, shape3, shape4) %>% 
  mcmc_dens()

post %>% select(lm1) %>% mcmc_trace()

library(coda)

coda::raftery.diag(posterior_samples(post))

> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     4        8092  3746          2.16     
>  sill    6        8524  3746          2.28     
>  range   6        8566  3746          2.29     
>  smooth  <NA>     <NA>  3746            NA     
>  shape1  3        4095  3746          1.09     
>  shape2  3        4095  3746          1.09     
>  shape3  4        8144  3746          2.17     
>  shape4  3        4163  3746          1.11     
>  shape5  3        4520  3746          1.21     
>  shape6  3        4095  3746          1.09     
>  shape7  3        4029  3746          1.08     
>  shape8  6        8532  3746          2.28     
>  shape9  3        4095  3746          1.09     
>  shape10 6        8342  3746          2.23     
>  shape11 4        6654  3746          1.78     
>  shape12 3        4267  3746          1.14     
>  shape13 3        4061  3746          1.08     
>  shape14 3        4267  3746          1.14     
>  shape15 3        4061  3746          1.08     
>  shape16 2        3897  3746          1.04     
>  shape17 2        3961  3746          1.06     
>  shape18 6        11313 3746          3.02     
>  shape19 3        4061  3746          1.08     
>  shape20 4        8092  3746          2.16     
>  shape21 3        4061  3746          1.08     
>  shape22 3        4028  3746          1.08     
>  shape23 3        4095  3746          1.09     
>  shape24 3        4267  3746          1.14     
>  shape25 3        4157  3746          1.11     
>  shape26 2        3865  3746          1.03     
>  shape27 2        3929  3746          1.05     
>  shape28 6        7792  3746          2.08     
>  shape29 3        4095  3746          1.09     
>  shape30 2        3994  3746          1.07     
>  shape31 3        4410  3746          1.18     
>  shape32 3        4028  3746          1.08     
>  shape33 3        4095  3746          1.09     
>  shape34 4        8298  3746          2.22     
>  shape35 6        8736  3746          2.33     
>  shape36 2        3994  3746          1.07     
>  shape37 3        4128  3746          1.10     
>  shape38 3        4128  3746          1.10     
>  shape39 3        4061  3746          1.08     
>  shape40 3        4095  3746          1.09     
>  chain   69075    69075 3746         18.40

# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.shape[,c(-4,-9)]), 
#                             as.mcmc(mc2$chain.shape[,c(-4,-9)]), 
#                             as.mcmc(mc3$chain.shape[,c(-4,-9)])))
# 
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.shape[,1]), 
#                             as.mcmc(mc2$chain.shape[,1]), 
#                             as.mcmc(mc3$chain.shape[,1])))

BMu1.mcmc<-mcmc.list(as.mcmc(mc1$chain.shape[,c(-4,-9)]), 
                     as.mcmc(mc2$chain.shape[,c(-4,-9)]),
                     as.mcmc(mc3$chain.shape[,c(-4,-9)]))

summary(BMu1.mcmc)

> 
> Iterations = 1:10000
> Thinning interval = 1 
> Number of chains = 3 
> Sample size per chain = 10000 
> 
> 1. Empirical mean and standard deviation for each variable,
>    plus standard error of the mean:
> 
>           Mean      SD  Naive SE Time-series SE
> lm1     0.7418 0.42890 0.0024763      0.0045351
> sill    0.1937 0.32464 0.0018743      0.0050592
> range   3.5989 5.38595 0.0310958      0.0908870
> shape1  0.6459 0.09745 0.0005626      0.0007165
> shape2  0.2402 0.10805 0.0006238      0.0009879
> shape3  0.6073 0.09190 0.0005306      0.0006901
> shape4  0.6922 0.09533 0.0005504      0.0007864
> shape6  0.6417 0.10743 0.0006202      0.0008532
> shape7  0.4316 0.08297 0.0004790      0.0005468
> shape8  0.5020 0.09171 0.0005295      0.0006737
> shape9  0.4534 0.10082 0.0005821      0.0007489
> shape10 0.5147 0.10265 0.0005926      0.0007973
> shape11 0.6614 0.09924 0.0005730      0.0007437
> shape12 0.8658 0.13977 0.0008070      0.0012274
> shape13 0.4066 0.10210 0.0005895      0.0007181
> shape14 0.4597 0.10748 0.0006205      0.0007699
> shape15 0.5024 0.10674 0.0006163      0.0007605
> shape16 0.4802 0.10774 0.0006221      0.0007219
> shape17 0.5337 0.11206 0.0006470      0.0007879
> shape18 0.6569 0.11091 0.0006403      0.0008731
> shape19 0.6658 0.09729 0.0005617      0.0007213
> shape20 0.9772 0.11750 0.0006784      0.0012123
> shape21 0.5361 0.10946 0.0006320      0.0008032
> shape22 0.6036 0.11916 0.0006880      0.0009378
> shape23 0.5628 0.12417 0.0007169      0.0010802
> shape24 0.5525 0.10463 0.0006041      0.0007478
> shape25 0.7334 0.10058 0.0005807      0.0008044
> shape26 0.5391 0.10400 0.0006004      0.0007936
> shape27 0.6670 0.09224 0.0005325      0.0007783
> shape28 0.8440 0.10508 0.0006067      0.0009787
> shape29 0.3711 0.10235 0.0005909      0.0007183
> shape30 0.5511 0.11441 0.0006605      0.0008423
> shape31 0.5225 0.10638 0.0006142      0.0008025
> shape32 0.4321 0.06945 0.0004010      0.0005309
> shape33 0.7575 0.09803 0.0005660      0.0007805
> shape34 0.4100 0.08273 0.0004776      0.0005993
> shape35 0.9596 0.11721 0.0006767      0.0011816
> shape36 0.6524 0.11323 0.0006537      0.0008618
> shape37 0.6166 0.11355 0.0006556      0.0008263
> shape38 0.5961 0.09975 0.0005759      0.0007665
> shape39 0.3269 0.11645 0.0006723      0.0009777
> shape40 0.8475 0.11604 0.0006700      0.0009448
> chain   2.0000 0.81651 0.0047141      0.0000000
> 
> 2. Quantiles for each variable:
> 
>            2.5%     25%     50%    75%   97.5%
> lm1     0.21272 0.54460 0.64613 0.8139  1.8975
> sill    0.02498 0.05159 0.08874 0.1944  1.0693
> range   0.25977 0.76101 1.64603 4.0021 19.5252
> shape1  0.46191 0.57955 0.64278 0.7097  0.8461
> shape2  0.02828 0.16409 0.24133 0.3159  0.4457
> shape3  0.43191 0.54578 0.60458 0.6658  0.7940
> shape4  0.51797 0.62601 0.68734 0.7506  0.8940
> shape6  0.44422 0.56619 0.63709 0.7135  0.8661
> shape7  0.27640 0.37399 0.42936 0.4868  0.6012
> shape8  0.32423 0.44009 0.50070 0.5637  0.6837
> shape9  0.25684 0.38609 0.45436 0.5197  0.6524
> shape10 0.32920 0.44343 0.50876 0.5784  0.7334
> shape11 0.46621 0.59598 0.66084 0.7268  0.8572
> shape12 0.61774 0.76700 0.85804 0.9559  1.1608
> shape13 0.21081 0.33785 0.40406 0.4735  0.6119
> shape14 0.25519 0.38571 0.46042 0.5305  0.6767
> shape15 0.30317 0.42997 0.49838 0.5696  0.7258
> shape16 0.28151 0.40469 0.47570 0.5500  0.7029
> shape17 0.32500 0.45801 0.53085 0.6045  0.7658
> shape18 0.45875 0.57873 0.65075 0.7275  0.8896
> shape19 0.47342 0.60086 0.66421 0.7317  0.8617
> shape20 0.76663 0.89403 0.97158 1.0527  1.2233
> shape21 0.33952 0.45975 0.52967 0.6072  0.7665
> shape22 0.37798 0.52251 0.60138 0.6814  0.8406
> shape23 0.32039 0.47840 0.56323 0.6486  0.8061
> shape24 0.34594 0.48262 0.55170 0.6218  0.7588
> shape25 0.53804 0.66647 0.73268 0.7997  0.9355
> shape26 0.34906 0.46804 0.53500 0.6029  0.7594
> shape27 0.49903 0.60412 0.66181 0.7250  0.8598
> shape28 0.64710 0.77226 0.84098 0.9125  1.0580
> shape29 0.17279 0.30149 0.37058 0.4387  0.5759
> shape30 0.33797 0.47217 0.54881 0.6280  0.7807
> shape31 0.31713 0.45309 0.52170 0.5908  0.7373
> shape32 0.29837 0.38406 0.43187 0.4788  0.5692
> shape33 0.57353 0.69037 0.75296 0.8216  0.9599
> shape34 0.25091 0.35412 0.40783 0.4646  0.5771
> shape35 0.74463 0.87792 0.95493 1.0362  1.2030
> shape36 0.44284 0.57410 0.64716 0.7265  0.8855
> shape37 0.40489 0.53909 0.61311 0.6907  0.8530
> shape38 0.41228 0.52667 0.59225 0.6608  0.8029
> shape39 0.10069 0.24825 0.32669 0.4045  0.5541
> shape40 0.63246 0.76711 0.84226 0.9239  1.0875
> chain   1.00000 1.00000 2.00000 3.0000  3.0000

xyplot(BMu1.mcmc)

densityplot(BMu1.mcmc)                             #Densidades

plot(BMu1.mcmc)

# layout(matrix(1:12, 3,4));  

traceplot(BMu1.mcmc)

#Grafica de autocorrelacion
autocorr.plot(BMu1.mcmc, auto.layout = TRUE, ask =F)

coda::acfplot(BMu1.mcmc)

 # gelman.diag(BMu1.mcmc)
# 
# gelman.plot(BMu1.mcmc)

geweke.diag(BMu1.mcmc)

> [[1]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>       lm1      sill     range    shape1    shape2    shape3    shape4    shape6 
> -0.541921 -0.911255 -0.635529  0.546572 -0.075900 -0.584468  1.271316 -0.874372 
>    shape7    shape8    shape9   shape10   shape11   shape12   shape13   shape14 
>  1.030963 -0.971523 -0.935714  0.256445  2.004252  0.695535 -0.413436  0.038511 
>   shape15   shape16   shape17   shape18   shape19   shape20   shape21   shape22 
>  1.790014  0.875726  0.739533  0.005974  1.728120  0.420641 -0.812484 -0.303818 
>   shape23   shape24   shape25   shape26   shape27   shape28   shape29   shape30 
> -1.669508 -0.062697  1.212304 -0.195292  0.512424  1.141296  0.221783  1.097566 
>   shape31   shape32   shape33   shape34   shape35   shape36   shape37   shape38 
>  0.015901  0.606067  2.901314  0.798967  0.277861 -0.220946 -0.381475  0.982637 
>   shape39   shape40     chain 
>  0.745817  0.855370       NaN 
> 
> 
> [[2]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>       lm1      sill     range    shape1    shape2    shape3    shape4    shape6 
> -0.541921 -0.911255 -0.635529  0.546572 -0.075900 -0.584468  1.271316 -0.874372 
>    shape7    shape8    shape9   shape10   shape11   shape12   shape13   shape14 
>  1.030963 -0.971523 -0.935714  0.256445  2.004252  0.695535 -0.413436  0.038511 
>   shape15   shape16   shape17   shape18   shape19   shape20   shape21   shape22 
>  1.790014  0.875726  0.739533  0.005974  1.728120  0.420641 -0.812484 -0.303818 
>   shape23   shape24   shape25   shape26   shape27   shape28   shape29   shape30 
> -1.669508 -0.062697  1.212304 -0.195292  0.512424  1.141296  0.221783  1.097566 
>   shape31   shape32   shape33   shape34   shape35   shape36   shape37   shape38 
>  0.015901  0.606067  2.901314  0.798967  0.277861 -0.220946 -0.381475  0.982637 
>   shape39   shape40     chain 
>  0.745817  0.855370       NaN 
> 
> 
> [[3]]
> 
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5 
> 
>       lm1      sill     range    shape1    shape2    shape3    shape4    shape6 
> -0.541921 -0.911255 -0.635529  0.546572 -0.075900 -0.584468  1.271316 -0.874372 
>    shape7    shape8    shape9   shape10   shape11   shape12   shape13   shape14 
>  1.030963 -0.971523 -0.935714  0.256445  2.004252  0.695535 -0.413436  0.038511 
>   shape15   shape16   shape17   shape18   shape19   shape20   shape21   shape22 
>  1.790014  0.875726  0.739533  0.005974  1.728120  0.420641 -0.812484 -0.303818 
>   shape23   shape24   shape25   shape26   shape27   shape28   shape29   shape30 
> -1.669508 -0.062697  1.212304 -0.195292  0.512424  1.141296  0.221783  1.097566 
>   shape31   shape32   shape33   shape34   shape35   shape36   shape37   shape38 
>  0.015901  0.606067  2.901314  0.798967  0.277861 -0.220946 -0.381475  0.982637 
>   shape39   shape40     chain 
>  0.745817  0.855370       NaN

#geweke.plot(BMu1.mcmc)
raftery.diag(BMu1.mcmc)

> [[1]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     3        4163  3746         1.11      
>  sill    4        5080  3746         1.36      
>  range   5        6077  3746         1.62      
>  shape1  3        4095  3746         1.09      
>  shape2  3        4095  3746         1.09      
>  shape3  3        4338  3746         1.16      
>  shape4  3        4163  3746         1.11      
>  shape6  3        4095  3746         1.09      
>  shape7  3        4031  3746         1.08      
>  shape8  2        3865  3746         1.03      
>  shape9  3        4095  3746         1.09      
>  shape10 2        3962  3746         1.06      
>  shape11 2        3962  3746         1.06      
>  shape12 3        4267  3746         1.14      
>  shape13 3        4061  3746         1.08      
>  shape14 3        4267  3746         1.14      
>  shape15 3        4061  3746         1.08      
>  shape16 2        3897  3746         1.04      
>  shape17 2        3962  3746         1.06      
>  shape18 2        3929  3746         1.05      
>  shape19 3        4061  3746         1.08      
>  shape20 4        4752  3746         1.27      
>  shape21 3        4061  3746         1.08      
>  shape22 3        4028  3746         1.08      
>  shape23 3        4095  3746         1.09      
>  shape24 3        4267  3746         1.14      
>  shape25 3        4147  3746         1.11      
>  shape26 2        3865  3746         1.03      
>  shape27 2        3929  3746         1.05      
>  shape28 3        4113  3746         1.10      
>  shape29 3        4095  3746         1.09      
>  shape30 2        3994  3746         1.07      
>  shape31 3        4410  3746         1.18      
>  shape32 3        4028  3746         1.08      
>  shape33 3        4095  3746         1.09      
>  shape34 3        4129  3746         1.10      
>  shape35 6        8738  3746         2.33      
>  shape36 2        3994  3746         1.07      
>  shape37 3        4129  3746         1.10      
>  shape38 3        4129  3746         1.10      
>  shape39 3        4061  3746         1.08      
>  shape40 3        4095  3746         1.09      
>  chain   <NA>     <NA>  3746           NA      
> 
> 
> [[2]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     3        4163  3746         1.11      
>  sill    4        5080  3746         1.36      
>  range   5        6077  3746         1.62      
>  shape1  3        4095  3746         1.09      
>  shape2  3        4095  3746         1.09      
>  shape3  3        4338  3746         1.16      
>  shape4  3        4163  3746         1.11      
>  shape6  3        4095  3746         1.09      
>  shape7  3        4031  3746         1.08      
>  shape8  2        3865  3746         1.03      
>  shape9  3        4095  3746         1.09      
>  shape10 2        3962  3746         1.06      
>  shape11 2        3962  3746         1.06      
>  shape12 3        4267  3746         1.14      
>  shape13 3        4061  3746         1.08      
>  shape14 3        4267  3746         1.14      
>  shape15 3        4061  3746         1.08      
>  shape16 2        3897  3746         1.04      
>  shape17 2        3962  3746         1.06      
>  shape18 2        3929  3746         1.05      
>  shape19 3        4061  3746         1.08      
>  shape20 4        4752  3746         1.27      
>  shape21 3        4061  3746         1.08      
>  shape22 3        4028  3746         1.08      
>  shape23 3        4095  3746         1.09      
>  shape24 3        4267  3746         1.14      
>  shape25 3        4147  3746         1.11      
>  shape26 2        3865  3746         1.03      
>  shape27 2        3929  3746         1.05      
>  shape28 3        4113  3746         1.10      
>  shape29 3        4095  3746         1.09      
>  shape30 2        3994  3746         1.07      
>  shape31 3        4410  3746         1.18      
>  shape32 3        4028  3746         1.08      
>  shape33 3        4095  3746         1.09      
>  shape34 3        4129  3746         1.10      
>  shape35 6        8738  3746         2.33      
>  shape36 2        3994  3746         1.07      
>  shape37 3        4129  3746         1.10      
>  shape38 3        4129  3746         1.10      
>  shape39 3        4061  3746         1.08      
>  shape40 3        4095  3746         1.09      
>  chain   <NA>     <NA>  3746           NA      
> 
> 
> [[3]]
> 
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95 
>                                                
>          Burn-in  Total Lower bound  Dependence
>          (M)      (N)   (Nmin)       factor (I)
>  lm1     3        4163  3746         1.11      
>  sill    4        5080  3746         1.36      
>  range   5        6077  3746         1.62      
>  shape1  3        4095  3746         1.09      
>  shape2  3        4095  3746         1.09      
>  shape3  3        4338  3746         1.16      
>  shape4  3        4163  3746         1.11      
>  shape6  3        4095  3746         1.09      
>  shape7  3        4031  3746         1.08      
>  shape8  2        3865  3746         1.03      
>  shape9  3        4095  3746         1.09      
>  shape10 2        3962  3746         1.06      
>  shape11 2        3962  3746         1.06      
>  shape12 3        4267  3746         1.14      
>  shape13 3        4061  3746         1.08      
>  shape14 3        4267  3746         1.14      
>  shape15 3        4061  3746         1.08      
>  shape16 2        3897  3746         1.04      
>  shape17 2        3962  3746         1.06      
>  shape18 2        3929  3746         1.05      
>  shape19 3        4061  3746         1.08      
>  shape20 4        4752  3746         1.27      
>  shape21 3        4061  3746         1.08      
>  shape22 3        4028  3746         1.08      
>  shape23 3        4095  3746         1.09      
>  shape24 3        4267  3746         1.14      
>  shape25 3        4147  3746         1.11      
>  shape26 2        3865  3746         1.03      
>  shape27 2        3929  3746         1.05      
>  shape28 3        4113  3746         1.10      
>  shape29 3        4095  3746         1.09      
>  shape30 2        3994  3746         1.07      
>  shape31 3        4410  3746         1.18      
>  shape32 3        4028  3746         1.08      
>  shape33 3        4095  3746         1.09      
>  shape34 3        4129  3746         1.10      
>  shape35 6        8738  3746         2.33      
>  shape36 2        3994  3746         1.07      
>  shape37 3        4129  3746         1.10      
>  shape38 3        4129  3746         1.10      
>  shape39 3        4061  3746         1.08      
>  shape40 3        4095  3746         1.09      
>  chain   <NA>     <NA>  3746           NA

# heidel.diag(BMu1.mcmc)

4 extra

bg <- bg %>% 
  dplyr::filter(Rain > 30)

bg %>% 
  plot_density(x = Rain,
               title = "(2018-2020) basic density plot general of cum_rain",
               fill = "salmon1")

table(bg$Station)
library(SpatialExtremes)
## Not run:
## Generate realizations from the model
n.site <- 4 # numero de sitios muestreados

n.obs  <- 11 # numero de observaciones tomada en cada sitio muestreado

coord  <- bg %>%                    # creamos la matrix de coordenadas
  select(Longitud, Latitud) %>%       # esta matrix de coordenadas solo incluye el
  unique()                            # numero de sitios donde se muestreo

coord  <- cbind(lon = coord$Longitud, lat = coord$Latitud) # matrix de coordenadas
# sin mediciones repetidas
# class(coord)
# coord <- cbind(lon = runif(n.site, -10, 10), lat = runif(n.site, -10 , 10))

# creamos un proceso gaussiano para cada parametro  de la GEV, esta parte se deja
# a elección del autor elegirlo.

gp.loc   <- SpatialExtremes::rgp(1, coord, "powexp", sill = 0.01, range = 10, smooth = 1)
gp.scale <- SpatialExtremes::rgp(1, coord, "powexp", sill = 0.001, range = 15, smooth = 1)
gp.shape <- SpatialExtremes::rgp(1, coord, "powexp", sill = 0.001, range = 20, smooth = 1)


locs   <- 26 + 0.5 * coord[,"lon"] + gp.loc
scales <- 10 + 0.2 * coord[,"lat"] + gp.scale
shapes <- 0.15 + gp.shape

SpatialExtremes::rgev(n.obs, locs[1], scales[1], shapes[1])

# Estimated Parameters:
#         xi         mu       beta 
#  0.3944641 34.1387472  4.7882714 

data <- matrix(NA, n.obs, n.site)
# 
# for (i in 1:n.site)
#   data[,i] <- SpatialExtremes::rgev(n.obs, locs[i], scales[i], shapes[i])
dim(bg)
table(bg$Station)

data[,1]  <-bg %>% 
  dplyr::select(Rain, Station) %>% 
  dplyr::filter(Station == "Arriaga") %>% 
  # slice(1:5) %>% 
  dplyr::select(Rain) %>% 
  dplyr::arrange(-Rain) %>% 
  slice(1:11) %>% 
  as.matrix()
data[,2]  <-bg %>% 
  dplyr::select(Rain, Station) %>% 
  dplyr::filter(Station == "Comitan") %>% 
  # slice(1:5) %>% 
  select(Rain) %>% 
  arrange(-Rain) %>% slice(1:11) %>% 
  as.matrix()
# data[,3]  <- bg %>% 
#   dplyr::select(cum_rain, Station) %>% 
#   dplyr::filter(Station == "scdlc") %>% 
#   # slice(1:5) %>% 
#   select(cum_rain) %>% 
#   arrange(-cum_rain) %>% # slice(1:61) %>% 
#   as.matrix()
data[,3]  <- bg %>% 
  dplyr::select(Rain, Station) %>% 
  dplyr::filter(Station == "Tapachula") %>% 
  # slice(1:5) %>% 
  select(Rain) %>% 
  arrange(-Rain) %>% slice(1:11) %>% 
  as.matrix()
data[,4]  <- bg %>% 
  dplyr::select(Rain, Station) %>% 
  dplyr::filter(Station == "Tuxtla") %>% 
  # slice(1:5) %>% 
  select(Rain) %>% 
  arrange(-Rain) %>%  slice(1:11) %>% 
  as.matrix()

# b2018[b2018$Station == "Arriaga", 3]
# data[,5] <- b2018[b2018$Station == "tuxtlag", 3]
# data[,1]  <- b2018[1:10,3]

# loc.form   <- y ~ lon
# scale.form <- y ~ lat
# shape.form <- y ~ 1

loc.form   <- y ~ 1
scale.form <- y ~ 1
shape.form <- y ~ 1

hyper <- list()
hyper$sills     <- list(loc = c(1,8), 
                        scale = c(1,1), 
                        shape = c(1,0.02))
hyper$ranges    <- list(loc = c(2,20), 
                        scale = c(1,5), 
                        shape = c(1, 10))
hyper$smooths   <- list(loc = c(1,1/3), 
                        scale = c(1,1/3), 
                        shape = c(1, 1/3))
hyper$betaMeans <- list(loc = 9, 
                        scale = 6, 
                        shape = 2)
hyper$betaIcov  <- list(loc =   solve(diag(c(10), 1, 1)),
                        scale = solve(diag(c(10), 1, 1)),
                        shape = solve(diag(c(0.13), 1, 1)))
## We will use an exponential covariance function so the jump sizes for
## the shape parameter of the covariance function are null.

prop <- list(gev = c(1.2, 0.08, 0.08), ranges = c(0.7, 0.8, 0.7), smooths = c(0,0,0))

start <- list(sills = c(4, .36, 0.009), 
              ranges = c(24, 17, 16), 
              smooths= c(1, 1, 1), 
              beta = list(loc = c(24), 
                          scale = c(7.31),
                          shape = c(0.54)))

# Estimated Parameters:
#   xi         mu       beta 
# 0.5304988 37.3221090  7.4043417 

mc1 <- latent(data, coord, 
             loc.form = loc.form, 
             scale.form = scale.form,
             shape.form = shape.form, 
             hyper = hyper, 
             prop = prop, 
             start = start,
             n = 10000, 
             burn.in = 5000, 
             thin = 15)
mc1

mc2 <- latent(data, coord, 
             loc.form = loc.form, 
             scale.form = scale.form,
             shape.form = shape.form, 
             hyper = hyper, 
             prop = prop, 
             start = start,
             n = 10000, 
             burn.in = 5000, 
             thin = 15)
mc2

mc3 <- latent(data, coord, 
             loc.form = loc.form, 
             scale.form = scale.form,
             shape.form = shape.form, 
             hyper = hyper, 
             prop = prop, 
             start = start,
             n = 10000, 
             burn.in = 5000, 
             thin = 15)
mc3
## End(Not run)

# a <- summary(mc)
# 
# a
# head(mc$chain.loc)
# head(mc$chain.scale)
# head(mc$chain.shape)
# head(mc$loc.dsgn.mat)
# head(mc$scale.dsgn.mat)
# head(mc$shape.dsgn.mat)

# tidybayes::add_residual_draws(mc$chain.loc[,6:10])

# mod1_sim <- coda::coda.samples(model = mc, 
#                          variable.names = chain.loc,
                         # n.iter = 5000)